##########################
# Replication file for the article:
# Title: "European Institutional Integration and the Educational Divide 
# in Support for the European Union"
# Authors: Sharon Baute & Tobias Tober
# Published: The British Journal of Political Science
# Date: April 2024
##########################

##########################
# Load packages
##########################

if (!require("brms")) install.packages("brms")
library(brms)
if (!require("arm")) install.packages("arm")
library(arm)
if (!require("tidybayes")) install.packages("tidybayes")
library(tidybayes)
if (!require("dplyr")) install.packages("dplyr")
library(dplyr)
if (!require("broom.mixed")) install.packages("broom.mixed")
library(broom.mixed)
if (!require("dotwhisker")) install.packages("dotwhisker")
library(dotwhisker)
if (!require("ggplot2")) install.packages("ggplot2")
library(ggplot2); theme_set(theme_bw())
if (!require("ggsurvey")) install.packages("ggsurvey")
library(ggsurvey)

##########################
# Load data
##########################

load("dat_macro.RData") # main dataset
load("dat_eurii.RData") # data for European institutional integration beginning in 1956

##########################
# Figure 1
##########################

dat_macro %>%
  dplyr::select(country, year, gmembrshp, wnation) %>%
  group_by(country, year) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1)) %>%
  ggplot(aes(year, perc)) + geom_line(size = 2) + 
  ylab("Percentage of support") + xlab("Year") + facet_wrap(~country) +
  ylim(0, 100) +
  theme(axis.title=element_text(size=22), 
        axis.text=element_text(size=16),
        strip.text.x = element_text(size = 18))

##########################
# Figure 2
##########################

dat_eurii %>%
  dplyr::select(country,country2,year,eurii) %>%
  ggplot(aes(year, eurii, group = country)) + geom_line(size = 2) + 
  ylab("European institutional integration") + xlab("Year") + 
  facet_wrap(~country) +
  theme(axis.title=element_text(size=22), 
        axis.text=element_text(size=16),
        strip.text.x = element_text(size = 18)) +
  geom_line(data = dat_eurii[,2:8], aes(year, eurii, group = country2), 
            color = "grey20", size = 0.9, alpha=.1)

##########################
# Figure 3
##########################

dat_macro %>%
  dplyr::select(eurii, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eurii, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1)) %>%
  ggplot(aes(eurii, perc)) + geom_point(size = 2) + geom_smooth(method = lm) +
  facet_wrap(~educats) +
  ylab("Percentage of support") + xlab("European institutional integration") +
  theme(axis.title=element_text(size=15), 
        axis.text=element_text(size=10),
        strip.text.x = element_text(size = 15))

##########################
# Figure 4
##########################

## Base model

# reduce dataset to relevant variables: no additional macro controls
vars <- c("gmembrshp","membrshp_ord","age","educats","sex","lrs","unmpl","married",
          "retr","eurii","country","eb","year","wnation")
dat1 <- dat_macro[vars]

# scale continuous variables
morph_vars <- c("age","lrs","eurii")
morph_vars_sc <- paste(morph_vars,"s",sep="_")
dat1 <- as.data.frame(na.omit(dat1))
for (i in 1:length(morph_vars)) {
  dat1[[morph_vars_sc[i]]] <- c(rescale(dat1[[morph_vars[i]]]))
}

# run model
m1 <- brm(gmembrshp ~ age_s + educats + sex + married
          + lrs_s + unmpl + retr + eurii_s + country
          +  (1|year) + (1|year:country),
          data = dat1, family = bernoulli,
          warmup = 500, iter = 1250, chains = 4, cores = 4,
          prior = c(set_prior("normal(0,1)", class = "b"),
                    set_prior("student_t(4,0,1)", class = "sd")))
coefs.m1 <- tidy(m1, effect = "fixed")
coefs.m1 <- coefs.m1[2:10,]
coefs.m1 <- coefs.m1[c(2:3,9,1,4:8),]

## Interaction model

# run model 
m2 <- update(m1, formula. = ~ . + educats:eurii_s,
             warmup = 500, iter = 1250, chains = 4, cores = 4,
             prior = c(set_prior("normal(0,1)", class = "b"),
                       set_prior("student_t(4,0,1)", class = "sd")))
summary(m2)
coefs.m2 <- tidy(m2, effect = "fixed")
coefs.m2 <- coefs.m2[c(2:10,25:26),]
coefs.m2 <- coefs.m2[c(2:3,9:11,1,4:8),]

## Model with country controls

# reduce dataset to relevant variables, including country controls
vars2 <- c("gmembrshp","age","educats","sex","lrs","unmpl","married",
           "retr","eurii","country","eb","year","unemp","sstran",
           "kof")
dat2 <- dat_macro[vars2]

# scale continuous variables
morph_vars2 <- c("age","lrs","eurii","unemp","sstran","kof")
morph_vars_sc2 <- paste(morph_vars2,"s",sep="_")
dat2 <- as.data.frame(na.omit(dat2))
for (i in 1:length(morph_vars2)) {
  dat2[[morph_vars_sc2[i]]] <- c(rescale(dat2[[morph_vars2[i]]]))
}

# run model 
m3 <- update(m2, formula. = ~ . + unemp_s + sstran_s + kof_s,
             newdata = dat2,
             warmup = 500, iter = 1250, chains = 4, cores = 4,
             prior = c(set_prior("normal(0,1)", class = "b"),
                       set_prior("student_t(4,0,1)", class = "sd")))
summary(m3)
coefs.m3 <- tidy(m3, effect = "fixed")
coefs.m3 <- coefs.m3[c(2:13,28:29),]
coefs.m3 <- coefs.m3[c(2:3,9,13:14,1,4:8,10:12),]

## Produce graph

# combine results
comb <- bind_rows(Controls=coefs.m3, Interaction=coefs.m2,
                  Base=coefs.m1, .id = "model") %>%
  filter(effect=="fixed")

# set fontface
fontface <- c("plain","plain","plain","plain","plain","plain",
              "plain","plain","plain","bold","bold","bold","bold",
              "bold","bold")

# plot
{dwplot(comb, dodge_size = 1, dot_args = list(aes(shape = model),
                                              size = 6)) +
    geom_vline(xintercept=0,lty=2) +
    xlim(-1,1) +
    scale_y_discrete(labels=c("age_s"="Age",
                              "sex1"="Female",
                              "married1"="Married",
                              "lrs_s"="Left-Right\nself-placement",
                              "unmpl1"="Unemployed",
                              "retr1"="Retired",
                              "eurii_s"="European inst. integration",
                              "educats1"="Medium education\n(Ref. Low)",
                              "educats2"="High education\n(Ref. Low)",
                              "educats1:eurii_s"="Medium education x\nEuropean inst. integration",
                              "educats2:eurii_s"="High education x\nEuropean inst. integration",
                              "unemp_s"="Unemployment rate",
                              "sstran_s"="Social transfers",
                              "kof_s"="Globalization")) +
    theme(axis.text = element_text(size=20),
          axis.text.y = element_text(face = fontface),
          legend.title = element_text(size=20),
          legend.text = element_text(size=20),
          legend.background = element_rect(colour="grey80")) +
    scale_color_grey(name = "Model:",
                     breaks=c("Base","Interaction",
                              "Controls"),
                     labels=c("Base"="Base",
                              "Interaction"="Interaction",
                              "Controls"="Country controls")) +
    scale_shape_manual(values=c(17,15,19),
                       name = "Model:",
                       breaks=c("Base","Interaction",
                                "Controls"),
                       labels=c("Base"="Base",
                                "Interaction"="Interaction",
                                "Controls"="Country controls"))} 

##########################
# Figure 5
##########################

# average marginal predicted probabilities of EU support conditional on
# education based on m1
tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$educats <- 0
sim_dat <- tmpdat
batch_size <- 5e4
sim_predictededu <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictededu[[length(sim_predictededu) + 1]] <- p
}
sim_predictededu <- data.frame(do.call(rbind, sim_predictededu))
plot_predictededu <- sim_predictededu %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$educats <- 1
sim_dat <- tmpdat
sim_predictededu1 <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictededu1[[length(sim_predictededu1) + 1]] <- p
}
sim_predictededu1 <- data.frame(do.call(rbind, sim_predictededu1))
plot_predictededu1 <- sim_predictededu1 %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$educats <- 2
sim_dat <- tmpdat
sim_predictededu2 <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictededu2[[length(sim_predictededu2) + 1]] <- p
}
sim_predictededu2 <- data.frame(do.call(rbind, sim_predictededu2))
plot_predictededu2 <- sim_predictededu2 %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))
eduplotdat <- data.frame(rbind(plot_predictededu,
                               plot_predictededu1, 
                               plot_predictededu2))
eduplotdat$var <- factor(0:2, levels = c(0,1,2),
                         labels = c("Low","Medium","High"))

eduplot <- ggplot(eduplotdat, aes(var, Estimate_name, group = 1)) +
  geom_errorbar(width = .2, aes(ymin = Q2.5_name, ymax = Q97.5_name)) +
  geom_point(shape=19, size=4) + ylim(0.48, 0.73) +
  xlab("Education") + ylab("Predicted probability of support") +
  theme(axis.title=element_text(size=22), 
        axis.text=element_text(size=16),
        strip.text.x = element_text(size = 18))

# average marginal predicted probabilities of EU support conditional on
# European institutional integration based on m1
tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
ivalues <- with(dat1, seq(from = min(eurii_s, na.rm = T),
                          to = max(eurii_s, na.rm = T),
                          length.out = 5))
isvalues <- with(dat1, seq(from = min(eurii, na.rm = T),
                           to = max(eurii, na.rm = T),
                           length.out = 5))
tmpdat$eurii_s <- ivalues[1]
sim_dat <- tmpdat
sim_predictedint <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictedint[[length(sim_predictedint) + 1]] <- p
}
sim_predictedint <- data.frame(do.call(rbind, sim_predictedint))
plot_predictedint <- sim_predictedint %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$eurii_s <- ivalues[2]
sim_dat <- tmpdat
sim_predictedint2 <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictedint2[[length(sim_predictedint2) + 1]] <- p
}
sim_predictedint2 <- data.frame(do.call(rbind, sim_predictedint2))
plot_predictedint2 <- sim_predictedint2 %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$eurii_s <- ivalues[3]
sim_dat <- tmpdat
sim_predictedint3 <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictedint3[[length(sim_predictedint3) + 1]] <- p
}
sim_predictedint3 <- data.frame(do.call(rbind, sim_predictedint3))
plot_predictedint3 <- sim_predictedint3 %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$eurii_s <- ivalues[4]
sim_dat <- tmpdat
sim_predictedint4 <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictedint4[[length(sim_predictedint4) + 1]] <- p
}
sim_predictedint4 <- data.frame(do.call(rbind, sim_predictedint4))
plot_predictedint4 <- sim_predictedint4 %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
tmpdat$eurii_s <- ivalues[5]
sim_dat <- tmpdat
sim_predictedint5 <- list()
for (i in seq(0, nrow(sim_dat), batch_size)){
  print(paste("# predicting from", i + 1, "to", i + batch_size))
  p <- fitted(m1, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
  sim_predictedint5[[length(sim_predictedint5) + 1]] <- p
}
sim_predictedint5 <- data.frame(do.call(rbind, sim_predictedint5))
plot_predictedint5 <- sim_predictedint5 %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))

intplotdat <- data.frame(rbind(plot_predictedint,
                               plot_predictedint2, 
                               plot_predictedint3,
                               plot_predictedint4,
                               plot_predictedint5))
intplotdat$var <- isvalues

intplot <- ggplot(intplotdat, aes(var, Estimate_name, group = 1)) +
  geom_line() +
  geom_ribbon(aes(ymin = Q2.5_name, ymax = Q97.5_name), alpha = 0.1) +
  xlab("European institutional integration") + ylab("") + ylim(0.48, 0.73) +
  theme(axis.title=element_text(size=22), 
        axis.text=element_text(size=16),
        strip.text.x = element_text(size = 18))

# combine both plots
eduint <- ggarrange(eduplot, intplot, labels = c("A", "B"),
                    font.label = list(size = 30))

##########################
# Figure 6
##########################

# average marginal predicted probabilities of EU support conditional on
# European institutional integration and education based on m2
tmpdat <- dat1[, c("age_s","educats","sex","lrs_s","unmpl",
                   "married","retr","eurii_s","country","year")]
jvalues <- with(dat1, seq(from = min(eurii_s, na.rm = T),
                          to = max(eurii_s, na.rm = T),
                          length.out = 10))
jsvalues <- with(dat1, seq(from = min(eurii, na.rm = T),
                           to = max(eurii, na.rm = T),
                           length.out = 10))

tmpdat$educats <- 0
batch_size <- 1e5
predprob1 = list()
for (j in jvalues) {
  print(paste("j =", j))
  sim_dat <- tmpdat
  sim_dat$eurii_s <- j
  sim_predicted <- list()
  for (i in seq(0, nrow(sim_dat), batch_size)){
    print(paste("# predicting from", i + 1, "to", i + batch_size))
    p <- fitted(m2, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
    sim_predicted[[length(sim_predicted) + 1]] <- p
  }
  sim_predicted <- data.frame(do.call(rbind, sim_predicted))
  sim_predicted$j <- j
  predprob1[[length(predprob1) + 1]] <- sim_predicted
}
edu0 <- do.call(rbind, predprob1)
edu0dat <- edu0 %>%
  group_by(j) %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))
edu0dat <- as.data.frame(edu0dat)

tmpdat$educats <- 1
batch_size <- 1e5
predprob1 = list()
for (j in jvalues) {
  print(paste("j =", j))
  sim_dat <- tmpdat
  sim_dat$eurii_s <- j
  sim_predicted <- list()
  for (i in seq(0, nrow(sim_dat), batch_size)){
    print(paste("# predicting from", i + 1, "to", i + batch_size))
    p <- fitted(m2, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
    sim_predicted[[length(sim_predicted) + 1]] <- p
  }
  sim_predicted <- data.frame(do.call(rbind, sim_predicted))
  sim_predicted$j <- j
  predprob1[[length(predprob1) + 1]] <- sim_predicted
}
edu1 <- do.call(rbind, predprob1)
edu1dat <- edu1 %>%
  group_by(j) %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))
edu1dat <- as.data.frame(edu1dat)

tmpdat$educats <- 2
batch_size <- 1e5
predprob1 = list()
for (j in jvalues) {
  print(paste("j =", j))
  sim_dat <- tmpdat
  sim_dat$eurii_s <- j
  sim_predicted <- list()
  for (i in seq(0, nrow(sim_dat), batch_size)){
    print(paste("# predicting from", i + 1, "to", i + batch_size))
    p <- fitted(m2, ndraws = 1000, newdata = sim_dat[(i + 1):min(nrow(sim_dat), (i + batch_size)), ])
    sim_predicted[[length(sim_predicted) + 1]] <- p
  }
  sim_predicted <- data.frame(do.call(rbind, sim_predicted))
  sim_predicted$j <- j
  predprob1[[length(predprob1) + 1]] <- sim_predicted
}
edu2 <- do.call(rbind, predprob1)
edu2dat <- edu2 %>%
  group_by(j) %>%
  summarize_at(vars(Estimate, Est.Error, Q2.5, Q97.5), 
               list(name = mean))
edu2dat <- as.data.frame(edu2dat)
edu2dat$js <- jsvalues
edu1dat$js <- jsvalues
edu0dat$js <- jsvalues

predprob <- ggplot(edu2dat, aes(js, Estimate_name)) +
  geom_line() + 
  geom_ribbon(aes(ymin = Q2.5_name, ymax = Q97.5_name), alpha = 0.1) +
  geom_line(data = edu1dat, aes(js, Estimate_name)) +
  geom_ribbon(data = edu1dat, aes(ymin = Q2.5_name, ymax = Q97.5_name), 
              alpha = 0.1) +
  geom_line(data = edu0dat, aes(js, Estimate_name)) +
  geom_ribbon(data = edu0dat, aes(ymin = Q2.5_name, ymax = Q97.5_name), 
              alpha = 0.1) +
  ylab("Predicted probability of support") +
  xlab("European institutional integration") +
  annotate("text", x = 75, y = 0.73, label = "High", size = 6) +
  annotate("text", x = 75, y = 0.625, label = "Medium", size = 6) +
  annotate("text", x = 75, y = 0.545, label = "Low", size = 6) +
  theme(axis.title=element_text(size=22), 
        axis.text=element_text(size=16),
        strip.text.x = element_text(size = 18))

##########################
# Figure 7
##########################

# Italy 
italy <- dat_macro[dat_macro$country=="Italy",]
italypre <- subset(italy, eb>290 & eb<350 & eb!=311 & eb!=321 & eb!=341)
italytreat <- subset(italy, eb>349 & eb<391 & eb!=351 & eb!=371 & eb!=381)

italy_92 <- italytreat %>%
  dplyr::select(eb, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eb, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1))

italypremean <- italypre %>%
  dplyr::select(eb, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eb, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1)) %>%
  group_by(educats) %>%
  summarize(meanpre = mean(perc))

it1 <- italytreat %>%
  dplyr::select(eb, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eb, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1)) %>%
  {
    ggplot(., aes(eb, perc, group = educats)) + 
      geom_line(aes(linetype = educats), size = 1.2) +
      scale_linetype_manual("Education", 
                            labels = c("Low", "Medium", "High"),
                            values=c("dashed", "dotted", "solid")) +
      scale_x_continuous(labels = .$eb, breaks = .$eb) +
      ylab("Percentage of support") + xlab("Wave") +
      geom_vline(aes(xintercept = 370), linetype = "dashed") +
      theme(axis.title=element_text(size=15), 
            axis.text=element_text(size=15),
            strip.text.x = element_text(size = 13),
            legend.title = element_text(size=15),
            legend.text = element_text(size=13),
            legend.key.width=unit(1.5,"cm"))
  }

italyplot <- italy_92 %>%
  full_join(italypremean, by = "educats") %>%
  mutate(as_mean = perc - meanpre) %>%
  {
    ggplot(., aes(eb, as_mean, group = educats)) + 
      geom_line(aes(linetype = educats), size = 1.2) +
      scale_linetype_manual("Education", 
                            labels = c("Low", "Medium", "High"),
                            values=c("dashed", "dotted", "solid")) +
      scale_x_continuous(labels = .$eb, breaks = .$eb) +
      ylab("Abnormal change in support") + xlab("Wave") +
      geom_vline(aes(xintercept = 370), linetype = "dashed") +
      theme(axis.title=element_text(size=15), 
            axis.text=element_text(size=15),
            strip.text.x = element_text(size = 13),
            legend.title = element_text(size=15),
            legend.text = element_text(size=13),
            legend.key.width=unit(1.5,"cm"))
  }

# UK
uk <- dat_macro[dat_macro$country=="United Kingdom",]
ukpre <- subset(uk, eb>290 & eb<350 & eb!=311 & eb!=321 & eb!=341)
uktreat <- subset(uk, eb>349 & eb<391 & eb!=341 & eb!=351 & eb!=371 & eb!=381 & eb!=391 & eb!=411)

uk_92 <- uktreat %>%
  dplyr::select(eb, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eb, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1))

ukpremean <- ukpre %>%
  dplyr::select(eb, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eb, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1)) %>%
  group_by(educats) %>%
  summarize(meanpre = mean(perc))

uk1 <- uktreat %>%
  dplyr::select(eb, gmembrshp, wnation, educats) %>%
  drop_na(educats) %>%
  group_by(eb, educats) %>%
  summarize(perc = 100*(weighted.mean(as.numeric(gmembrshp), wnation, na.rm=T)-1)) %>%
  {
    ggplot(., aes(eb, perc, group = educats)) + 
      geom_line(aes(linetype = educats), size = 1.2) +
      scale_linetype_manual("Education", 
                            labels = c("Low", "Medium", "High"),
                            values=c("dashed", "dotted", "solid")) +
      ylab("Percentage of support") + xlab("Wave") +
      geom_vline(aes(xintercept = 370), linetype = "dashed") +
      theme(axis.title=element_text(size=15), 
            axis.text=element_text(size=15),
            strip.text.x = element_text(size = 13),
            legend.title = element_text(size=15),
            legend.text = element_text(size=13),
            legend.key.width=unit(1.5,"cm"))
  }

ukplot <- uk_92 %>%
  full_join(ukpremean, by = "educats") %>%
  mutate(as_mean = perc - meanpre) %>%
  {
    ggplot(., aes(eb, as_mean, group = educats)) + 
      geom_line(aes(linetype = educats), size = 1.2) +
      scale_linetype_manual("Education", 
                            labels = c("Low", "Medium", "High"),
                            values=c("dashed", "dotted", "solid")) +
      scale_x_continuous(labels = .$eb, breaks = .$eb) +
      ylab("Abnormal change in support") + xlab("Wave") +
      geom_vline(aes(xintercept = 370), linetype = "dashed") +
      theme(axis.title=element_text(size=15), 
            axis.text=element_text(size=15),
            strip.text.x = element_text(size = 13),
            legend.title = element_text(size=15),
            legend.text = element_text(size=13),
            legend.key.width=unit(1.5,"cm"))
  }

# combine four plots
estudy <- ggarrange(it1, uk1, italyplot, ukplot, labels = c("IT","UK","IT","UK"),
                    font.label = list(size = 25), label.x = -0.02)